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ABSTRACT 

Upcoming surveys for galaxy clusters using the Sunyaev-Zel'dovich effect are poten- 
tially sensitive enough to create a peculiar velocity catalog. The statistics of these 
peculiar velocities are sensitive to cosmological parameters. We develop a method to 
explore parameter space using N-body simulations in order to quantify dark matter 
halo velocity statistics which will be useful for cluster peculiar velocity observations. 
We show that mass selection bias from a kinetic Sunyaev-Zel'dovich velocity catalog 
forecasts rms peculiar velocities with a much more complicated f2 m dependency than 
suggested by linear perturbation theory. In addition, we show that both two-point 
functions for velocities disagree with linear theory predictions out to ~ 40 ft. -1 Mpc 
separations. A pedagogical appendix is included developing linear theory notation with 
respect to the two-point peculiar velocities functions. 

Key words: cosmology: theory - cosmology: observation - cosmology: large-scale 
structure of the Universe - galaxies: clusters: general 



1 INTRODUCTION 



The growth of galaxy and galaxy cluster peculiar ve- 
locities provides information on the growth of struc- 
ture in the gravitational instability paradigm. Over the 
last decade, cosmic velocity fields were an active area 
of research, b oth in observation and in t heory (e.g., 
iBahcall fc Obi <|l99fih : IStrauss fc Willickl <|l995h l. Bulk pe- 
culiar velocities of galaxies and clusters were modelled as 
tracers of the background dark matter velocity field and 
were frequently used to constrain cosmolog ical parameters 
even f airly recently (^B jndlg_etalJ, |200lj; iFgidrnanejL^lJ, 
| 2003l: lJuszkiewicz et all Il999t ISheth fc Diaferiol 1200 3: 
iPeel fc Knoxll2002ft . 

As measured by the number of relevant papers and con- 
ferences, there has been a slightly falling interest in velocity 
work, in part due to observational limitations. The funda- 
mental plane method of measuring galaxy peculiar veloci- 
ties, based on methods such as Tully-Fisher and D n -cr, is 
limited by relati ve intrinsic errors, w hich grow as a percent- 
age of distance Jjacobv et all Il992ft . For a highly selected 
subsample, the errors can be as low as ~ 10 per cent, but in 
general the error is closer to 15-20 per cent of the distance. 
This has limited the direct use of velocities for cosmology 



to redshifts of 



0.024 A rou 



70 h' 1 Mpc (e.g- lBridle et alJ lioOll) ') (throughout this pa- 



a comoving distance of 
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per, h is the dimensionless hubble parameter such that to- 
day, H = h 100 kms" 1 Mpc" 1 ). 

In contrast, peculiar velocities derived from doppler- 
shifted ('kinetic') Sunyaev-Zel'dovich effect (kSZ) spectra 
are subject to entirely different systematic and intrinsic er- 
rors. The kSZ effect probes the hot gas within the cluster 
and represent s a noisy estimate of the cluster's bulk motion 
( Holdeil l2004lk Complex motions of the hot gas (cold fronts, 
coolin g flows, etc.) further increase the noise llNaeai et all 
12003ft . In addition, measuring the kSZ signal, which only re- 
flects the motion of the innermost p art of intrac l uster gas, 
is a challenging observational effort llKnox et al 1 l2004t ?); 
some recent experiments have achieved limited success (?). 
Nevertheless, as the kSZ effect is not redshift limited in prin- 
ciple (errors may grow indirectly with the distance as a func- 
tion of cluster-parent halo evolution), it is a very promising 
technique for constraining parameters. 

With the view that kSZ observational difficulties can 
be mitigated by sheer numbers and clever signal separation 
techniques, the question remains as to whether we are ap- 
plying the right theoretical velocity models. Typically, lin- 
earized first order perturbation theory ("linear theory" in 
this paper) has been used for large scale velocity fields, al- 
though modes at the galaxy scale (~ Mpc) are non-linear 
and even inter-galaxy scales (~10 Mpc) are quasilinear for 
2 ~ 0. This has been justified by only relying on the theory 
at larger scales where it is reasonable to invoke the stable 
clustering regime. 

Yet clusters are assumed to exhibit large fluctuations 



2 A. C. Peel 



at large scales and are clearly biased samples in the linear 
regime. How does this affect the model? In other words, what 
does selection bias do to the statistics of these velocities? In 
light of the growing body of work on non-li near halo evo- 
lution in simulations us i ng th e 'halo model' llMo fc WhiteL 
Il996l : ISheth fc Torment I1999T) . peculiar velocities are due 
for a similar detailed examination. This is especially true 
for so-called 'precision cosmology' efforts when constraining 
parameters such as fi m and as. This might seem obvious in 
dealing with galaxy peculiar velocities. But it is also true for 
the streaming motions of galaxy clusters, whose rarity im- 
plies that environment dependence, biasing, and selection 
effects may be more important than for galaxies. 

For instance, simulated clusters (large mass dark matter 
halos) s how rms pecul i ar vel ocities that depart from linear 
theory dColberg et all feoOOt) . An earlier paper found that 
linear theory predictions were 'somewhat lower' than N- 
body results and noted the disagreeme nt between simulated 
and l inear theory two-point functions JCroft & Efstathiou, 
Il995l) . Reasonable attempts to explain the excess rms of 
peculiar velocities using simulations have bee n published 
jsheth fc Diaferiol l200lt> : iHamana et al.l <2003l) '). although 
not specifically as functions of cosmological parameters. In 
this paper, we examine the full two-point velocity functions 
through N-body simulations while varying Q m ■ Understand- 
ing the behaviour of these functions will be necessary for 
observations to yield constraints on parameters. 

In §2, we provide a brief pedagogical discussion of lin- 
ear theory peculiar velocities and discuss why linear theory 
is overly simplistic. We address the peak-background split 
approach and examine how selecting over the peaks in the 
density field affects theoretical predictions. In §3 we discuss 
our simulations. In §4 we summarize our results. We dis- 
cuss our results in §5. Our conclusions in §6 also include a 
discussion of the usefulness of our approach for parameter 
forecasting in general. An Appendix is included in this pa- 
per which outlines how to calculate the two-point correlation 
functions for velocities of peaks. 



2 LINEAR THEORY 

2.1 The two-point correlation tensor 

We der i ve the two-point velocity functions in real space (e.g., 
iGorskil il988h ). (For mathematical elaboration, see also the 
Appendix in this paper.) This is crucial to properly con- 
struct the velocity correlation matrix used in any constraint 
analyses, as cross-correlations between velocities must be 
taken into account. 

By "linear theory", we mean that we begin with an 
initially Gaussian distributed field and evolve it in the lin- 
ear regime of the gravitational instability paradigm via 
perturbation theory within a standard Friedman- Walker- 
Robertson universe. 

With that simplification, the continuity equation 

k ■ v k = iS k (1) 

implies that the curl-free v k will grow as the time derivative 
of the density field, where 8k is the comoving mode of the 
density contrast 8p/p and v k is the Fourier velocity com- 
ponent parallel to that mode. The overdot is the conformal 
time derivative. 



z=0, R = 5 Mpc/h 
= 0.35 
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Figure 1. Linear theory velocity correlations of a Gaussian dis- 
tributed field smoothed with a 5h~ 1 Mpc window function — all 
other parameters are held fixed for a flat ACDM universe. From 
top to bottom: 9\\ , cj>n ($11 is discussed in §2.3); solid, dotted 
and dashed lines show f2 m =0.35, 0.3 and 0.25 respectively. 



The two-point statistic encompasses correlations be- 
tween two vectors, so we have the nine-element tensor: 

V ij (r) = (v i (x)v j (x + r)) = 1,2,3) (2) 

where we invoke isotropy and homogeneity so that ^ can 
only depend on the comoving distance r = \r | between pairs 
of velocities at comoving positions x — r\ and x + r — V2- 
The () brackets refer to an ensemble average. It is straight- 
forward to derive the two-point radial velocity function. On 
average, there are only two nontrivial correlations: one for 
the components of the velocities parallel to the line between 
them (^n), and one for the components perpendicular to 
that line ($±). We obtain: 

*(r) = + (*n - * ± )ff (3) 

where I is the identity tensor and r is the unit vector along 
r. 

Fig-Qshows how <Pm and depend on comoving dis- 
tance between two points for three different flat ACDM cos- 
mologies (parameters other than Q m fixed; $11 is discussed 
below in §2.3). Note the dependence on Q m is fairly degen- 
erate with normalization such as by 0$. In addition, these 
functions have been convolved with smoothing tophat win- 
dow function to account for the velocity of an extended sec- 
tion of the field, rather than for a point particle. 

For a given pair of radially projected velocities separated 
by angle 9 on the sky, v r ( , yi,ri), Vr(72, ^2), the two-point 
correlation is: 

*12 = (VrlVra) = 71 ■ * ■ 72 

= ¥j.coB0+(¥||-*j.)/(0,ri,ra) (4) 

where: 
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cos a = 71 • 72 
and 

f(0,ri,r 2 ) 



(rf + r 2 )cos6 - rir 2 (l + cos 2 9) 
r\ + r\ — 2r\T2 cos 6 



(•») 



(6) 



For the extreme cases of either two positions lined up along 
the line of sight: 



and ri ^ r-2 f 1 

#12 = #11 



(7) 



or for two positions at the same radial distance, but sepa- 
rated by a small angle: 

8 but small and r\ — f — > — — 

#12= (l-£)*±-£*||«*L. (8) 

The zero-lag for either two-point function is related to 
the rms velocity: 

a 2 v = <« 2 > = (v(x) ■ v(x)} = Ei{^ 2 > = 3# ( || or _L)(0) (9) 

where we note that the literature often uses a 2 and {v 2 ) 
interchangeably. 



2.2 Arguments against linear theory 

The naive idea that linear theory in the field will predict 
the velocities of galaxies should be suspect, though this 
has been the norm in the p ast (with some notable excep- 
tions, e.g.. lMa fc Frvl i2002l) . ?). Galaxies are not only non- 
linear objects themselves, their velocities are often respond- 
ing strongly in the non-linear regime since perhaps as many 
as ~ 20 per cent of them are in bound groups. Further- 
more, they most likely represent biased objects compared 
to the general behaviour of the dark matter background. 
Data from large scale surveys has been consistent with mod- 
elling thi s bias as approximately linear for L* (and dimmer) 
galaxies ( Seli ak fc WarrenL l2004h . though simulations sug- 
gest the bias is exp ected to increase for larger sized haloes 
jSeliak et aTL 12004) . 

Clusters are on average inherently rare objects as mod- 
elled by any Press-Schecter type formalism. Although a 
10 14 ft -1 Mq cluster began in a comoving volume of radius 
R ~ 7/i -1 Mpc (ACDM), few such volumes have clusters. 
In fact, a cluster is found at late times (on average) in a ra- 
dius ~4 times that size, i.e., a volume 64 times the original 
source volume for the cluster. This implies that galaxy clus- 
ters as a selected sample should have a bulk motion which is 
responding to long wavelength modes which have not under- 
gone collapse and are therefore well-modelled in the linear 
approximation. 

The statistical rarity of high mass haloes is presumed 
to be a source of a suppression factor to their rms ve- 
locity as predicted b y the excursion hierarchy approach of 
iBardeen et alJ dl98rJ) . Modeling clusters as originating from 
the 3-cr (or greater) end of the density peak distribution, the 
peak rms velocity is given in linear theory by: 



r 2 p (R) 



where: 



tJtp { k,t)k 2 "W 2 (kR) 



(11) 



= D 2 ( V ) 



k 2 dk 
2tt 2 



\5 k , \ 2 k 2n W 2 (kR). 



aZ(R) 



(10) 



Implicit in the final form above is the linear theory assump- 
tion that S(x,rj) is separable; thus 8{x,rf) = D(rj)8o(x) and 
we choose the growth function D to be normalized to one 
today. (D = a for an Q m — 1 cosmology). This assumption 
is reasonable if we smooth on large enough scales to ignore 
nonlinear processing on small scales, which is the purpose of 
the window function W(kR). The window function is usu- 
ally chosen to be the Fourier transform of either a Gaussian 
or tophat envelope. When the time (t or 77) is not speci- 
fied, D, or P(k) assume values at redshift zero, i.e., today. 
For examples of this notation, the cosmological parameter 
o"8 = o"o(8 h~ x Mpc), and the <j v of Eqs. © and 110H above 
is oc <T_i. 

The rms velocity of rare, massive objects is apparently 
suppressed by a few to t en percent compared to background 
as implied by Eq. lUTHl . Colberg et alJ J200CT) show that al- 
though this statement agrees for the velocities of the den- 
sity peaks at early times, large mass haloes evolved to low 
redshift which formed around these peaks actually have a 
higher rms velocity than linear theory would have predict 
by as much as 40 per cent. This discrepancy is observed 
in our simulations even at moderate redshifts (z Z 0.6) in 
cluster formation and will be discussed in §5 below. See the 
Appendix for more detailed calculations to predict peak- 
peak velocity correlations as a natural extension to Eq. 1101 : 
however, it is fairly clear from simulations that the peak- 
background split approach is not accurate at late times, as 
evidenced for example with the rms peak velocity above. 

How does one take into account the bias for dense 
objects, like galaxies and clusters, caused by their being 
in an overdense region? In general, an object's peculiar 
velocity is greatly affected by its environment: haloes in 
overdense re gions typically move faster than those i n less 
dense ones llColberg et all l200d : ISheth fc DiaferioL l200ll : 
lHamana et alll2003fl . Although the average intracluster dis- 
tance is large, the likelihood of finding a cluster near another 
is high, and of finding a cluster near a large overdensity of 
galaxies and groups is very high. If environment plays an 
important role for the evolution of galaxy peculiar veloci- 
ties, it must therefore play an even more important one for 
clusters. 



2.3 Momentum correlations 

A heuristic way to see the effect of this selection bias is to 
re-examine Eq. (0 within the context of linear theory. The 
full form, after separating out the background solution and 
regardless of the amplitude of 5(x) is: 

S + V- [(l + 5)v] =0 . (12) 
This suggests we might consider the statistic: 

# = (v(x)v(x + r)) 

= {v(x){l + 8{x))v(x + r)(l+5(x + r)) (13) 
= (V1V2) + (viv 2 5i) + (V1V2S2) + (V1S1V2S2) 
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where the subscripts *1' and '2' correspond to the arguments 
x and x + r respectively. This is a 'momentum correlation', 
i.e., weighting the velocity by the density in the region. 

Evaluating Eq. (I13P using arguments within linear the- 
ory solely to build intuition (i.e., three-point functions are 
zero, ignoring the connected part of the four-point function, 
etc.) leads to: 

*(r)=*_L(l+eW)i+((*||-*x)(l + ?W) + *||)^(14) 
with $x an d $11 as before, and: 
$|| = {viS2){v 2 Si) 

= - ° lD ^ lD2 dk |4,o| 2 fcji(fcr)^ 2 (fcJ?)|(15) 

£(r) is the usual Fourier transform of the power spectrum: 

y=-\5 k , \ 2 jo(kr)W 2 (kR). (16) 

Fig-Qshows the behaviour of |$| at z~0. 

So how does Eq. I14H compare to the unweighted model, 
Eq. iJ^J? The extra factor of £ boosts the correlations for ve- 
locities at short (~ 10 h^ 1 Mpc) separations. The new factor 
of $n which is zero at both zero lag and large separations 
boosts the anticorrelations in for a characteristic separa- 
tion. This is not presumed to be exact, but will aid discussion 
of results in §5. 

By definition, linear theory becomes invalid when 5 ~ 1 
which forces the equations in Fourier space to mix modes. 
This is where numerical simulations become essential. 



2.4 Previous efforts regarding the velocity rms 

Extensive w o rk b y ISheth fc Diaferid (1200 ll) and 
lHamana et alJ <l2003f) comparing halo velocities in simula- 
tions to linear theory showed that the local environment of 
a halo had a heavy influence on the evolution of its velocity. 
Specifically, the bias predicted by the halo model suggested 
that halo velocities would likely be boosted as if they had 
evolved in a higher- Q m universe (see Fig. because haloes 
are typically found in overdense regions. 

In one approach iSheth fe Diaferid. 2001), it was sug- 
gested that a typical halo speed today would be related to 
the linear growth velocity (evolved from redshift 20 until 
today) boosted by the local density in the region: 

vo = (i + sr^^l V20 a?) 

D(yo) 

where S is smoothed over a region of radius R using a Gaus- 
sian window function. They found that fi was naturally tied 
to the choice in smoothing radius, and for their simulations, 
fit: 

^- a6 ^10?M P c) - < 18 > 

Following this guideline, the second approach 
llHamana et all 120031) found a similar result phrased 
as the rms velocity (essentially the same statistic: see the 
final paragraph in §2.1 above): 

°Lio(M, S) = [l + 5{R local )] 2 » {Rl °«> l) al{M) (19) 



although neither group found the velocities to have a strong 
dependence on halo mass. 

In the first work, the choice of 10 h~ Mpc as a reference 
scale for smoothing is motivated by the fact that it roughly 
represents the transition from linear to non-linear regimes 
today as measured by u(R) ~ 1. According to the second 
work, deciding how to choose Ri oca i is primarily an ansatz. 

The velocity statistics in previous papers investigated a 
wider range of dark matter halo masses. We will focus on the 
statistics of only the largest mass haloes by imposing a mass 
cutoff sugg ested from model s of the Sunyaev-Zel'dovich ef- 
fect iCarlstrom et all 120021) . In addition, we will examine 
this selection effect on the full two-point functions, rather 
than simply the zero- lag. 



3 SIMULATIONS 

3.1 ART code 

We used an Ada ptive Refinement Tree (ART) N-body code 
teravtsovt I1999T) which, like its predecessor Particle Mesh 
(PM) N-body codes, integrates trajectories of collisionless 
particles by solving the Poisson equation. Unlike PM codes, 
ART allows for a hierarchy of refinement meshes where col- 
lapsed objects require more resolution. 

ART employs standard particle-mesh techniques to 
compute acceleration grids in order to advance particle co- 
ordinates and velocities in time. A regular cubic grid covers 
the entire computational volume and defines the initial min- 
imum resolution of the simulation. This grid is then refined 
where the density contrast is higher to form higher resolu- 
tion sub-meshes in those regions of interest. The main com- 
putational loop of the integration consists of: (1) density 
assignment for all existing meshes; (2) running the gravi- 
tational solver; (3) routine updating particle positions and 
velocities; (4) modifications to the mesh hierarchy. 

3.2 Halo finder algorithms 

The ART codes we used produce files of particle positions 
and velocities which were subsequently analysed for the pres- 
ence of haloes. The basic problem of halo finding in a simu- 
lation is that there are no clear boundaries for haloes. There 
is no single perfect algorithmic definition of a group or mass 
of a group. 

Many halo finding algorithms exist, but tend to fall into 
two categories: the friends of friends type (FoF) linked-list 
type approaches where particles are identified with a halo 
if they are within a cer tain chosen distance of each other 
fefstath iou et alT Il985l) ; and overdensity methods such as 
DENMAX which calculates the density as a function of a 
grid and identifies to which local m aximum each particle 
belongs feertschinger fc Gelbl 1199 ll). We used a rel atively 
recent method named HOP lEisCTafcein fc Hutill998t) which 
follows the logic behind overdensity methods yet includes 
'hopping' to nearest neighbors a la FoF methods. Instead 
of calculating the density on a grid, a density is associated 
with each particle. Then a search is conducted for the highest 
density nearest neighbor until a particle is its own densest 
neighbor. All particles which trace to the same such parti- 
cle are grouped. A followup 'regrouping' then reunites any 
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sufficiently bound haloes which happen to contain two (or 
more) local maxima such that the initial hopping misiden- 
tified them as separate haloes. 

3.3 Virial radius 

After HOP was used to find the haloes, a crude spherical 
overdensity method was applied to restrict the statistics to 
different cutoff radii. Real measurements of cluster peculiar 
velocities via the kSZ effect will be restricted to the baryons 
at the core but are likely to represent at best t he bulk motion 
of pa rticles 'trapped' within the virial radius llHolder et al.L 
2001). We define the virial radius by beginning nearest the 
central overdensity of a halo and including particles at every 
increasing radii until the overdensity within that radius is 
180 times the background density. 



where a 2 (R) is the usual smoothing of the power spectrum 
with a window function W(kR), and the number of haloes 
of mass M at a redshift z is: 

Jl(M, z)dM = 2^vf(v(M))dp (22) 
where 

v = S c /a(M) (23) 

and 5c is the critical value of a spherical overdensity at 
turnaround time. The universality referred to is due to the 
functional form of / and is not as useful for our purposes as 
Eq. J23 above. For 300 high mass clusters (3x 10 14 h" 1 M Q ) 
at a redshift of z ~ 0.6, we required a fairly large volume 
of (850 /i -1 Mpc) 3 . The steepness of the halo mass function 
would then guarantee ~ 10 3 clusters with mass greater than 
2 x 10 14 h' 1 M B at z ~ 0.6. 



3.4 Preliminaries 

The first questions to answer using N-body simulations were 
to determine: (1) how many high-mass haloes (presumably 
hosting clusters) were needed; (2) how big the simulated 
volume should be; and (3) how much mass resolution was 
required for each halo. This phase was completed using ap- 
proximately 10,000 hours of processor time on the the IBM 
SP computer, 'Seaborg', at the National Energy Resource 
Computing Facility at Lawrence Berkeley National Labora- 
tory. 



3. 4-1 Number of high mass haloes (clusters) 

We relied on linear theory predictions as a rough guide in 
determining the number of high mass (M Z 3 x 10 14 M@) 
haloes we would need to achieve an error variance on the 
order of a percent for Q m . 

For effectively uncorrelated cluster velocities, from, e.g., 
a very sparse survey, we estimate the expected error variance 
on fi m . From a measurement of iV clusters with their pecu- 
liar velocity variance represented by the zero-lag value of 
either two-point function (^o(zi)): 



(A!!,; 



E 



\ dQ.m 



2(*o(z 4 )+< 1 



(20) 



where the last equality assumes TV clu sters with a: 
*o all at 2=1, and dln^ /dn m ~ 5 JPeel fc Knoxl 12002 
Thus, on the order of 1000 clusters would be apparently suf- 
ficient to constrain fi m to a few percent. Current and future 
cluster surveys expect to detect on the or der of 10,000 clus- 
ters t hrough the Sunyaev-Zel'dovich effect dCarlstrom et al.L 
2002). 



3.4-3 Number of particles per halo 

We simulated a volume with the same initial conditions but 
with three different mass resolutions to see how many parti- 
cles were needed to resolve halo velocities. For this conver- 
gence test, we used smaller boxes of 150 Mpc per side 
for speed. At this size, we expected very few haloes above 
10 14 h' 1 M Q at a redshift of z ~ 0.6. 

Beginning with 256 3 cells, we used number of particles 
64 3 , 128 3 and 256 3 and tracked the velocities of the top five 
haloes as they became more resolved, as well as the rms of 
the entire population of haloes. From this convergence test, 
it became clear that the number of dark matter particles 
required to resolve a velocity was approximately 70, which 
meant that for an 850 Mpc sized box, (256) 3 particles 
would be sufficient to resolve the velocities of the largest 
haloes, i.e., those haloes most likely to contain massive clus- 
ters. Eight 425 h~ l Mpc per side simulations would cover the 
same volume and run much faster, but with the loss of the 
k = 27r/850 /i Mpc -1 mode. We ran multiple 425 h' 1 Mpc 
per side simulations with 128 3 particles and compared re- 
sults with one 850 h' 1 Mpc using 256 3 particles. We found 
that losing the low- A; mode (fc = 27r/850ft Mpc -1 ) had a 
negligible effect on velocity statistics. 

For our chosen resolution of 128 3 particles realized in 
ten (425 h' 1 Mpc) 3 volume boxes, m p = 3.0 x 10 12 h' 1 M 
(ACDM, Q m =0.3). This means that any halo identified with 
approximately 60 particles (1.8 x 10 14 h^ 1 M©) is a halo ca- 
pable of hosting a cluster. The simulations were run using 
three values of Q m (0.25, 0.3, 0.35) in flat ACDM cosomolo- 
gies with all other parameters fixed on the UK National 
Cosmology Supercomputer in Cambridge. Each realization 
took less than one week on eight Altix 3700 Itanium2 proces- 
sors. (Eight Altix processors was approximately ideal for this 
number of particles and timesteps.) We discuss the broad 
usefulness of this approach in §6. 



3-4-2 Volume 

To determine the necessary volume t o find these haloes, we 
followed Jenkins (I Jenkins et all 1200 If) fitting formula for the 
'universal mass function': 

f(M) = 0.315 exp(-|lna _1 +0.61j 3 ' 8 ) (21) 



4 RESULTS 

Table 1 summarizes the the rms peculiar velocity and 
number of haloes found above the cutoff mass as a func- 
tion of cosmology in our total volume (10 x(425) 3 ~ 
(915 h~ x Mpc) 3 ). For brevity, we will initially display 
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Table 1. Number and Peculiar Velocities of High Mass Haloes 
from N-body Simulations 
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0.25 


0.3 


0.35 


0.25 


0.3 


0.35 


0.0 


9344 


12618 


16541 


532 


539 


534 


0.25 


6211 


8235 


10530 


554 


552 


543 


0.667 


2320 


3009 


3790 


555 


542 


536 


1.5 


109 


131 


146 


563 


522 


462 



our results in the context of only one cosmology, a flat 
ACDM with f2 m =0.3. We first consider the Gaussianity of 
the one-dimensional velocity distribution (and the related 
Maxwellian speed distribution) to justify error estimates. 
Then we show the zero lag portion of the two-point func- 
tions: the rms peculiar velocities as a function of mass and 
density. Following this, we reveal the primary result of this 
work: the redshift evolution of the two-point functions for 
haloes above a cutoff mass. At that point, we will also in- 
troduce our results from two other cosmologies, flat ACDM 
with f2 m =0.25 and 0.35. We conclude this section by show- 
ing the bulk velocity history of the particles which make up 
the zero redshift haloes above a cutoff mass. We discuss and 
explain the results in §5. 

Only haloes with a minimum mass of 1.8 x 10 14 /i _1 
Mq (60 particles in our fiducial fl m —0.3 runs) were used to 
ensure reliable velocities (see §3.4.3 above). Our conclusions 
will not depend on this low mass estimate for a cluster capa- 
ble of producing a measurable Sunyaev-Zel'dovich effect; in 
fact, they would be more dramatic for higher mass cutoffs. 
After identifying haloes by their member particles as de- 
scribed above, we calculated the bulk flow peculiar velocity 
by averaging particle velocities within the virial radius, de- 
fined as the spherical radius for which the halo is 180 times 
the background density. The resulting velocity vector is then 
used to calculate all halo velocity statistics. 

To understand how appropriate linear theory is for cer- 
tain scales (and as a consistency check on our simulations), 
we also track the average velocities and two-point functions 
for the field. To do this, each simulation is partitioned into 
equal-sized boxes and the particle velocities within each box 
are averaged to create a 'bulk' peculiar velocity. The process 
is repeated for different smoothing lengths. We coin these 
partitions 'miniboxes.' The statistics of these minibox bulk 
velocities represented a simulated version of linear theory, 
convolved with cubic window functions. 

This partitioning was repeated for different nominal lin- 
ear sizes of L — 8-128 h~ Mpc by powers of two. Since the 
simulations were of linear size 425 /t" 1 Mpc, these lengths 
were rescaled due to roundoff to 8.02, 16.3, 32.7, 70.8, and 
142 /i" 1 Mpc so that the entire volume would be sampled. 
These values correspond to tophat-sphere smoothing vol- 
umes with radii of 4.98, 10.1, 20.3, 43.9 and 87.9 h' 1 Mpc. 
This was tested by examining the theory calculation with 
either a boxcar window function or a tophat window func- 
tion. In either case the same theoretical results were ob- 
tained when the volumes of the respective cubes or spheres 
in real space were equivelant. Each minibox was then as- 
signed a density based on the number of particles found 
inside: Si — rii/n — 1. Attempts to look at miniboxes of 



smaller extent were limited by the spatial resolution of the 
simulation. 

Three things should be mentioned about how the lin- 
ear theory was calculated. First, the apparent dependence 
on mass shown in the downward curvature of the linear the- 
ory is due to convolving the integral in Fourier space with 
a tophat window function W(kR), with a mass associated 
with the comoving volume 4nR 3 p/3 = M in contras t to th e 
smaller smoothing radii used bv lCroft fc Efstathioul il995t) . 
Secondly, the suppression factor mentioned in §2.2 above is 
completely ignored as we are comparing linear theory only 
to the entire particle field through the minibox statistics. 
Finally, the integration limits in Fourier space were cho- 
sen to match the simulation: the lowest frequency associ- 
ated with the size of the box represents the lower bound on 
the fc-integrals and the upper limit value was chosen to re- 
flect the Nyquist frequency. Initially, the Nyquist frequency 
was simply 2tt/425 x 128/2 h Mpc" 1 « 1 h Mpc" 1 . While 
the dynamic mesh was resolved by as many as six times at 
late redshifts, those new distances represent collapsed re- 
gions and do not force much of an increase in the Nyquist 
frequency for integrating the linear power spectrum. Since 
we integrate with a window function, the effect of using 1 h 
Mpc" 1 vs. (2 6 ) = 64 h Mpc" 1 for the upper limit is negli- 
gible. 

4.1 Gaussianity of velocity distributions 

We examine velocity distributions at redshifts z = for 
haloes and miniboxes in order to understand relevant con- 
fidence intervals for v 2 . If the set of one dimensional veloc- 
ities {vi ■ x,Vi ■ y,Vi ■ z}, is Gaussian, then {v 2 } should 
fit a Maxwellian distribution. Figs. [5] and [3] below show 
that Gaussian (and related Maxwellian) distributions work 
reasonably well for both field (as represented by mini- 
boxes) and halo one-dimensional velocities (and related 3-d 
speeds), even for differ ent mass ranges. The bias seen in 
ISheth fc Diaferiol teOOll) between smaller and larger haloes 
is much less apparent for masses above ~ 10 14 h' 1 M Q . Both 
distributions are fit by a Maxwellian with one-dime nsion al 
velocity dispersion of 311 kms" 1 , i.e., an rms of {v 2 } = 
V3a v i — 539 kms" 1 . The shot noise from rarity at higher 
redshift for a large Sunyaev-Zel'dovich survey clearly im- 
plies a transition to a Poisson distribution when the data is 
binned in redshift. This is also true from our work, but is 
not shown here for brevity. 

There is a high velocity tail, as predicted by 
ISheth fc Diaferiol I^OOlf l. However, for our large number of 
cluster-sized haloes, the majority of the distribution is still 
well enough fitted to a Maxwellian to justify the l-cr errors 
in the figure below given by: 

2 1 ■-' (24) 



{V ) 



In the case of the two-point function, this becomes: 

(*i,iiW> = ^l*x,nW*x,n(o)l ■ 



(25) 



4.2 Velocity rms 

Fig. [I] shows the rms velocities of miniboxes and haloes 
vs. linear theory predictions as functions of mass at differ- 
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Figure 2. One-dimensional velocity distributions (left) 
([vx, Vy, v z ]) and speed distribution (right) for miniboxes 
smoothed over 8.02 Mpc (top) and 32.7 h^ 1 Mpc at redshift 
0. Matching Gaussian and Maxwellian distributions for the 
one-dimensional velocities and speed respectively are shown as 
dotted lines. Recall that v rms = \fZa v \. 
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Figure 3. As above, but for haloes binned in two mass ranges at 
redshift 0. 



ent redshifts. In the mass ranee of interest (10 14 - fewxlO 15 ) 
linear theory and the minibox velocity statistics are in excel- 
lent agreement. Although not shown for brevity, the minibox 
two-point functions were also in agreement with linear the- 
ory. 

This figure also shows the discrepancy between rms ve- 
locities of linear theory and haloes. This is in rough agree- 



Figure 4. Peculiar velocity rms as a function of mass: boxes show 
halo statistics; triangles show minibox statistics; solid lines show 
linear theory. The circle shows the average halo statistic for the 
entire mass range of interest with insignificant statistical 1-sigma 
error bar. 



ment with the result that linear theory underpredicts the 
rms of cluster peculiar velocities as represented by massive 
dark matter haloes by a large percent. For example, averaged 
across the mass range, the predicted value of halo rms pecu- 
liar velocity is 539±3 kms -1 at redshift zero, whereas linear 
theory predicts an average of 415 km s~ i n that range, a dis- 
crepan cy o f 30 per cent. In agree ment with lSheth fc Diaferiol 
( 2001) and Ha mana et alJ a2003r) . we find essentially no mass 
dependence for the rms peculiar velocities (within l-cr). 
Since linear theory does predict a weak dependence on mass 
(larger objects should be slower as shown by the convex line 
in Fig.[IJ the discrepancy found bv lColberg et all ll2000l) ac- 
tually worsens for the largest mass haloes, although clearly 
this is of less statistical significance. 



4.3 Two-point functions 

We show the two-point function redshift dependence it 
two ways. First, motivated by obervations, we examine the 
haloes above a cutoff mass which have collapsed at redshifts 
z=0.0, 0.25 and 0.667, shown in Figs. |S| |7| respectively. 
We also examine the velocity history of the particles which 
will be within the virial radius of z=0 haloes above the mass 
cutoff. We do this by calculating particles' 'bulk' velocity as 
if they had formed haloes already. We only use those par- 
ticles which will be within the virial radius at z=0. Those 
results at redshifts z=0.25, 0.667 and 19 are shown in Figs.|H| 
l9l 1101 respectively. 

For Figs. through [5] we use linear theory smoothed 
over R = 10 h' 1 Mpc as a benchmark for comparison, 
though we do not a priori expect halo velocity statistics 
to be in agreement with linear theory. This raises the ques- 
tion as to what the appropriate comparison smoothing scale 
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Figure 5. Halo velocity two-point correlations (perpendicular 
and parallel components) at redshift z = 0.0 for haloes above the 
mass cutoff (see text). The dotted, short-dashed and long-dashed 
lines show linear theory and the triangles, squares and pentagons 
show simulated haloes for f2 m = 0.25, 0.3, and 0.35 respectively. 
Note the high value for simulations at the zero lag (where ty± 
must equal ^ii). 



should be. The average distance at z—19 of the particle far- 
thest from the eventual halo centre was calculated to be 
11.81, 11.23 and 10.64 ft" 1 Mpc for cosmologies fi m =0.25, 
0.3, and 0.35 respectively. Fig. llOl shows these values rather 
than R = 10/i _1 Mpc which would reflect larger values of 
the two-point function at zero lag (by approximately 5 per 
cent). (Larger values of the smoothing scale suppress the 
linear two-point functions more at short-distances.) 



4-3.1 as functions of redshift and cosmology 

Examining the two-point functions for different £l m values 
at different redshifts, four crucial results stand out. First, 
looking at Table 1, it is essentially impossible to discrimi- 
nate between these three different cosmologies at zero lag at 
redshift z=0, in contrast to the behaviour predicted by linear 
theory JPeel fc Knoxl f2002l. Second, at higher redshifts, the 
rms peculiar velocity dependence on Q m for these three val- 
ues is exactly opposite to that expected from linear theory 
(sec Fig.0. Third, the behaviour of the parallel component 
shows heavy influence from infall for r < 30 ft -1 Mpc at any 
redshift for which these massive haloes exist. In particular, 
the extreme anticorrelation seen in Fig.[7|at r ~ 4 ft" 1 Mpc 
is comparable (or larger) in magnitude to the zero lag value 
shown. Fourth, for the three values of Q m the perpendicu- 
lar components alone seem to reflect linear theory closely in 
behaviour if not in amplitude. 



Figure 6. as above, but for z = 0.25. 
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Figure 7. as above, but for z = 0.667. 



4-3.2 ^ as functions of history 

Tracking the history of the particles which will assemble the 
largest haloes by z~0 (as well as those at higher redshift), 
it is clear that the parallel component shown in Figs.|H|and 
does not reflect the infall as strongly as seen in Figs. © 
and[7| Note also that the parallel component shows remark- 
ably little difference in behaviour (apart from amplitude) 
between z=19 and 2=0.667. Prior to the assembly of the 
largest haloes, the gross behaviour of their particles is nearly 
fixed to track linear theory. Only just as and after a halo viri- 
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Figure 8. Velocity two-point correlations at z = 0.25 (perpendic- 
ular and parallel components) for the groups of particles which 
will be in the haloes shown in the Fig. [F] The dotted, short- 
dashed and long-dashed lines show linear theory and the trian- 
gles, squares and pentagons show simulated haloes for Q m = 0.25, 
0.3, and 0.35 respectively. Note the high value for simulations at 
the zero lag (where \I>j_ must equal Wii). 



alizes does the strong infall out to r < 30/i _1 Mpc become 
apparent. This is discussed in §5 as a coincidence in timing. 

In Fig. I1UI the rms peculiar velocity is smaller than 
linear theory predictions (and would be even smaller if we 
had retained R = 10 /i" 1 Mpc smoothing). This suppression 
is comp arable to the excursio n heirarchy suppression pro- 
posed bv lBardeen et all dl986l) mentioned in §2.2 above. The 
jump between the perpendicular component at r Z 2h _1 
Mpc and the rms value at r=0 is just as abrupt as it is at 
later redshifts. Some of this may be explained as an artefact 
of resolution (recall the average interparticle separation is 
3 ft -1 Mpc). Nonetheless, the perpendicular component for 
the particles destined for large haloes is actually less like 
linear theory here then at any subsequent time. 



5 DISCUSSION 

There are at least three different factors which help explain 
the differences between cluster-sized dark matter halo pe- 
culiar velocity two-point functions at different values of Q m 
and their differences from the two-point functions predicted 
by linear theory. Specifically, there are two selection biases 
to consider as well as the effect of dark energy domination. 
The first selection bias (also considered in previous work) 
is that regions which harbor the seeds of large dark matter 
haloes are by definition overdense. This helps explain why 
the rms peculiar velocity is 30 per cent higher than predicted 
by linear theory at any redshift in which one would find such 
a virialized large dark matter halo. 

The second, more subtle selection bias strongly affects 



Figure 9. as above, but for z = 0.667. 
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the constraining power of a galaxy cluster-based set of veloc- 
ity observations. The rarity of large haloes is sensitive to a 
combination of f2 m and as . If one fixes the value of as but de- 
creases the value of Q m , there is an effective transfer of power 
from small scales to large ones. This creates deeper initial 
potential wells at the (now rarer) largest scales and thereby 
increases the acceleration of halo velocities. The largest mass 
haloes will preferentially be found near these large rare fluc- 
tuations. If we could calculate the velocity statistics using all 
haloes down to galaxy size (many of which are not respond- 
ing to the extremely rare, deep, large-scale potentials), we 
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would find that there was a decrease in the velocity power 
at small scales and would recover the expected fl m depen- 
dence predicted by linear theory. However, the hot plasma 
necessary to produce a Sunyaev-Zel'dovich effect observa- 
tion requires a deep potential well and is therefore sensitive 
to a specific halo mass cutoff. This is why the rms peculiar 
velocities of the f2 m =0.25 haloes for a high mass cutoff were 
comparable or in excess of those for larger values of f2 m even 
at higher redshifts. The very slight mass dependence (which 
seems statistically insignificant in Fig.|lJ is exactly the issue 
when varying Q m . It is clear that if we were to use an ever 
larger mass cutoff in that figure, the average rms velocity 
would increase while the number of haloes used decreased. 

The anticorrelation in the parallel component is sim- 
ply a result of infall, a non-linear process to which the lin- 
ear theory is predictably b lind. This behaviour was noted in 
ICroft fc EfstathiouHl995D : in an effort to continue using lin- 
ear predictions, the authors adjusted the smoothing length 
in the linear model and selected velocities below a cutoff. 
They conclude this approach introduces a strong bias in pa- 
rameter determination. We examined this effect at different 
redshifts to determine how this discrepancy with linear the- 
ory has evolved. By redshift z=0.667, the few haloes massive 
enough to produce a Sunyaev-Zel'dovich effect are so rare as 
to have come from extremely overdense regions where col- 
lapse has been accelerated as if from a much higher Q. m 
universe. Consequently, infall is apparent out to large sepa- 
rations (~ 40/i _1 Mpc). However, at lower redshifts where 
the rarity of such regions decreases, the anticorrelation scale 
and its overall (absolute) magnitude also decrease because 
the simulation is beginning to enter the era of dark energy 
domination (simply represented in our work by a cosmolog- 
ical constant). The extreme anticorrelation seen in Fig. |7| 
gives way to a much milder anticorrelation only seen at the 
smallest separations by 2=0. Dark matter-cosmological con- 
stant equality occurs at 2=0.44, 0.33 and 0.23 for f2 m =0.25, 
0.3 and 0.35 respectively and marks the beginning of the end 
for new infall. Large modes which have not yet collapsed be- 
gin to decay during the acceleration phase, so volumes not 
undergoing gravitational collapse never will. 

Fig. demonstrates that the cosmological constant has 
become dominant more recently for the larger Q m values 
because the parallel anticorrelation at r ~ 2.5 ft -1 Mpc is 
still comparable (though negative) to about half the rms 
velocity. In contrast, for fi m =0.25, the growth of structure at 
these comoving scales has already halted and the expansion 
has accelerated sufficiently to begin a decay towards linear 
theory values, especially for r k, 5 fa -1 Mpc. 

Finally, we remark on the history of particles destined 
to be in large mass haloes by 2=0. The rms peculiar velocity 
derived from these sets of particles at 2=19 (Fig. HOI is lower 
than the linear theory prediction in agreement wi th th e 
peak-background split prediction of iBardeen et all |l986). 
In addition, the perpendicular and parallel components for 
separations r ^ 40 ft -1 Mpc are also lower in amplitude than 
linear theory predictions. However, as the haloes assemble 
in the heirarchical paradigm (Figs. |5] and |UJ, the rms pecu- 
liar velocity surpasses linear theory and the perpendicular 
components (from 2 ^ r ^ 40 h' 1 Mpc) evolve towards lin- 
ear theory values. The perpendicular component represents 
pairs of haloes responding only to some third large-scale fluc- 
tuation. It therefore probes the field more effectively than 



the parallel component in which the pairs' self-attraction at 
these scales dominates (on average) over gravitational at- 
tractions from other potential wells. Consequently, the per- 
pendicular component gradually recovers the behaviour pre- 
dicted by linear theory, modulo an overall boost in ampli- 
tude created by the mass selection bias mentioned above. 
In contrast, the parallel component gradually reflects the 
ever more common large scale infall until matter domina- 
tion ends. The effect is heuristically like the function $y 
mentioned in §2 above, though the scale occurs at smaller 
separation and the effect has a larger amplitude due to the 
non-linearity of the local density. 

By 2=0.25, when many of the 2=0 haloes have more 
than 60 particles within their virial radii, the parallel com- 
ponents begin to display the anticorrelation more appar- 
ent in Figs. I5JJ7] The coincidence in timing mentioned in 
§4.3.2 refers to the rarity of large scale potential wells com- 
pared to their amplitudes. As the largest mass haloes be- 
come more common, the infall felt by the ones still forming 
comes from more common smaller amplitude potential wells; 
consequently, there is less apparent infall in Fig.|5]than Fig.Q 
simply because the volume probed by the particles destined 
for haloes is much larger than the rare, high-amplitude fluc- 
tuation dominated volume probed by fully formed high mass 
haloes. 



6 CONCLUSIONS 

The purpose of this work was two-fold. The primary goal 
was to examine how fl m alone affects the velocity statistics 
of cluster-sized haloes using simulations as opposed to what 
linear theory predicts. The secondary purpose was to show 
that small (128 3 particle) N-body simulations are sufficient 
to characterize the statistics of cluster-sized dark matter 
haloes and furthermore are fast enough to allow an explo- 
ration of parameter space in a reasonable amount of time. 

We have shown that the dependence on fl m for the rms 
peculiar velocity for a fixed value of as and a fixed lower cut- 
off for the halo mass is counterintuitive to what linear theory 
predicts with a fixed erg. The peculiar velocities which de- 
velop for the largest mass haloes are faster for a smaller 
value of Q m while the universe is matter dominated because 
the largest modes must have a greater amplitude for a fixed 
value of as. However, since the decay of the (uncollapsed) 
largest modes occurs earlier for a smaller Q m , there is a 
convergence towards an indistinguishable rms peculiar ve- 
locity at 2 ~0 even though (in our work) f2 m varied by ~ 16 
per cent. Peculiar velocity catalogs from kinetic Sunyaev- 
Zel'dovich observations therefore have a much more com- 
plicated dependence on fi m than other data from, e.g., the 
Cosmic Microwave Background, or large-scale galaxy red- 
shift surveys. 

In addition, we have shown that the two-point veloc- 
ity functions (in real space) for cluster-sized dark matter 
haloes are as sensitive to these complicating biases as the 
rms peculiar velocity alone. The prospect of using the infor- 
mation from the two-point functions themselves to constrain 
fi m are currently hindered by two major observational fac- 
tors. First, the mass-based selection bias mentioned above 
affects the perpendicular component in a manner very sim- 
ilar to its effect on the rms velocity. Second, the sources of 
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noise mentioned in §1 are not guaranteed to be stochas- 
tic, yet are co mparable in amplitude to the rms signal 
(« 200 kms-^ ljKnox et, all Eool iNaeai et "all l200fl . 

The full cosmological sensitivity of a kinetic Sunyaev- 
Zel'dovich derived peculiar velocity catalog has not been 
explored. This is where our secondary goal will come into 
play. Rough scaling shows the potential for a reasonable 
search through parameter space. For example, based on clus- 
ter statistics using multiple small (128 3 particle) ART simu- 
lations on a current, fairly large supercomputer or multiple 
efficient clusters, we find: 



T 



(■A^ parameters \ 
5 J 



N, 



values / parameter 



-^realizations/ value \ / 8cpu/ realization 



10 



\ A^cpus 



(26) 



150 weeks. 



For instance, five parameters at three values each with ten 
realizations using only 160 CPUs would run in about two 
months. While the rest of Eq. 126H does scale linearly, the 
number of CPUs per realization is a non-linear function of 
the architecture of the computer and the performance of 
the simulation. The choice of 8cpu/reaiization is based on our 
simulations for this project and is only given as a guideline. 
Considering Moore's law, this an ever more conservative es- 
timate. This is not on par with the speed of a Monte-Carlo 
Markov Chain method. However, when simulations are the 
only reliably accurate way to explore parameter space, a 
simpler Fisher Matrix type approach is now conceivable for 
multiple parameters. 
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APPENDIX: TWO-POINT CORRELATIONS AT 
PEAKS 

This appendix is provided primarily as a pedagogical tool. 
The nonlinear evolution of clusters and the mass selection 
effect of kSZ observations renders the consideration of peak 
vs. background statistics moot here, especially since, as 
shown below and in past papers, the peak statistics predict 
suppressed velocities compared to background, i.e., general 
field velocities. However, future methods of observation and 



statistical analysis not mentioned in this paper may find the 
calculations presented below quite relevant. 

W e parallel the approach detailed in iBardeen et alJ 
( 1986) for calculating the peak statistics for a Gaussian dis- 
tributed linear density field, 8(x). The density gradients at 
the peak, rji(x p ) — ViS\ Xp , are zero (and nearly zero in 
some peak neighborhood) and the second derivatives of the 
density at the peak, C,ij(x p ) = ViV jd\ x , form a 3 x 3 sym- 
metric matrix with positive eigenvalues. The six indepen- 
dent entries are relabelled (a where {A — 1,2,3} corre- 
spond to the diagonal elements (originally £11,(22 and (33) 
and {A = 4, 5, 6} correspond to the off-diagonals (originally 
(23, Ci3 and (12). 

We shall adopt much of the notation found in 
IBardeen et alJ il986h . We replace (1,(2,(31 with x,y,z 
where: 



X = -((1 +(2 +( 3 )/o-2 

V = -(Ci - C 3 )/(2a 2 ) 



(27) 



-((l-2( 2 +( 3 )/(2t72) 



Then, 5 correlates only with x, not y, z, 77, nor (4,5,6- 
We also scale 8 by (To, and thus have constructed unitless 
variables with simple Gaussian widths: {v 2 = 8 2 /a 2 ) = 
(x 2 ) — 15{y 2 } = 5{z 2 ) — 1. In addition, if we ignore the 
velocities, we have a nearly symmetric auto-correlation ma- 
trix: the only non-zero, independent off-diagonal element is 
(vx) = o\l(olol). 

Of course, we are also concerned with the velocities. 
In linear perturbation theory, the divergence of the veloc- 
ity field is given by the time derivative of the density field 
(Eq. Q). Therefore, each component Vi only correlates with 
rji, and to estimate the peak rms velocity we only need to 
examine the 6x6 element matrix: 



M = 



(vv) {vrj) 
(rjv) (vv) 



&v 13 

Da 2 I 3 



Da 2 Is 
a 2 I 3 



(28) 



(I3 is simply the 3x3 identity matrix) . Inverting this matrix 
and looking at the quadratic form Q = d T M^ 1 d/2 which 
appears in the Gaussian pdf P oc e~® (d = [vi,rji]), while 
imposing the condition that rf(x p ) = leads to the expres- 
sion given m Eq. ab ove for the autocorrelations of peak 
velocities. 

Complicating matters, however, is the fact that 
many of the corresponding two-point correlations are 
not zero. One must reexamine the entire correla- 
tion matrix to see the effects on the peak veloc- 
ity two-point functions (or, in fact, any other two- 
point functions). The list of relevant variables is now 

d = [5i, Wi,^,^, 1/1, Zi,(ia,<5 2 ,V 2 ,T7 2 , £2,2/2, 22, (2A; A = 

4, 5, 6], where the subscripts {1, 2} refer to the points x and 
x + r, respectively. Our correlation matrix has now blos- 
somed to 26 x 26 elements 1 . 

If we continue to follow IBardeen et alJ lll986l) . the next 
step would be to rotate ( onto its principle axes. However, 
we cannot rotate both (1 and (2 simultaneously. Instead, it 
makes more sense to rotate the z axis to be in alignment 



1 "Even for the two-point function, the task of int egrating over all 
these variables is not pleasant to contemplate." IBardeen et all 
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with the line-of-sight between the two points, i.e., along r. 
This is what is done in evaluating Eq. ||3J. 

The full analytic forms for the linear perturbation the- 
ory velocity correlation functions are often given as: 



and: 
*7> 



Dit) 2 f k 2 dk |5 M ' 2 



(29) 



2tt 2 



k 2 



Mkr)+j 2 (kr)]W'(kR) 



= / dk\5 k , \ 2 [j (kr) - 2j 1 (kr)/kr]W 2 (kR) (30) 



E>l£>2 f k 2 dk \ s k,o\ 



3 -fe&^Mkr) 



2j 2 (kr)]W 2 (kR) 



where we have explicity added a —1 superscript: in general, 
both the parallel and perpendicular ^"s will be defined as 
above, with a k 2 " in the integral, much like the role of the 
n subscript for Eq. Hill . We also define: 

^\5k,o\ j m (kr) k n W 2 (kR). (31) 



C(r-) = DiDi 



Below, we will implicitly assume that all <f> and a func- 
tions are in terms of growth functions, not their time deriva- 
tives; thus we will use appropriate factors of /3 — D/D when- 
ever velocities appear. 

The second lines in Eqs. 1291 and 1301 make it clear that 
we can rewrite \P" above as Pi(3 2 (<f>o n + <^2 n )/3 and similarly 
= Pifh(4>l n - 20l")/3. Also note that £(r) = <$>{r) and 
a n = <Po n (0). In this way, all the auto- and two-point cor- 
relations are expressible as linear combinations of the set of 
functions 0J^(r). As mentioned above, we will always rotate 
the z-axis to be parallel to r. 

We solve for and list the various auto- and two-point 
correlations below, with subscript a — 1,2 for the two points 
x and x + r respectively. For clarity, we also label the 
Kronecker-delta function with a superscript K. For com- 
pleteness, the auto-correlations are repeated here: 



{VaVa) 


= 1 


{VaiVaj} 




(VaiTjaj) 


= S^p a a 2 /3 


(VaiVaj) 




{v a X a ) 


= di/(crocr 2 ) 


(XaXa) 


= 1; (VaVa) = 1/15; 


{C,aiC,aj) 


= Sfja 2 /15 

{i,3 = 4,5,6} 



all others being zero. For the following, we suppress the r 
argument as given, and rotate z || r. The two-point correla- 
tions are: 

(v2Vu) = Stepi^/ao (32) 

(viT]2i) = -<5i30l/°"O = ~{V2T}U) 

(vix 2 ) = 4>l/{o (j 2 ) 

{V1D2) = 4>l/(2(T()(T 2 ) 

(v\Zi) = —4>i/(2a a o 2 ) 
(vibi) = 0; {i = 4,5,6} 



(»i«a> = /3i/3 2 [*I 1 I+(*^ 1 -*I 1 )ff] 
(wii| a ) = ^ll+i^-^Drr] 

{V2V1) = /&[*51+ (33) 

(viiX 2 ) = -5 i3 P\4>\/(T2 

(vuV2) = 5f 3 pi(4>\/5- 3^/10) /<T2 
{vnz 2 ) = ~{vuy 2 } 

(wiiOy) = 0i(<A + <&)(S*83i + 65S? s )/5 {j = 4,5,6} 

(^1^2) = -&a<j>\/02 = —(TfriXi) (34) 

(r?iiJ/2) = *«(0i/5 - 3tfi/10)/ffa = -<i/w»i) 

{iliiz 2 ) = —{r)uy 2 ) = —(772^-81) 

faiiCw) = (0? + *s)(*£*£ + *£$)/5 0=4,5,6} 

(ansa) = </>o/cr| 

(K1J/2) = (j&Kloi) = ~(X1Z2) 

(y iy2 ) = [1901/140 - 0f/21 + <^/15]/<r 2 2 

< yi z 2 > = (0l/7 - 30|/28)/a 2 2 

(*iza> = [27^/140 + 2 /7 + <?f>o/5]/<7 2 2 

<Ci4<24> = [-4^/35 -^/21 + ^/15]/a 2 2 

= (C15C25) 

<Ci6C2 6 ) = [<^/35 + 2^/21 + <^/15]/a 2 2 

Note that most of the correlations are the same under r — > 
— r except those involving one T7 or one 1; which pick up 
a minus sign e.g., Eqs. 1)320 and 1)34)1 : and those involving 
only one v which forces D\ —+ D 2 and D 2 — > _Di, e.g., 
Eqs. 1321 and 1331 1. In the limit where we ignore correlations 
at distances large enough to reflect a significant difference 
between D\ and D 2 or between their time derivatives, this 
becomes a nitpicking concern. 

To understand peak velocity two-point correlations, we 
must now perform the following steps. First, we construct 
the required 26 x 26 matrix M using the above results and in- 
vert it. Next, we integrate the probability distribution given 
by P oc e' Q where Q = d T M~ 1 d/2 over all the variables 
except the velocities: over v\ and v 2 from the peak cutoff 
(such as ~ 3) to infinity; and over all values of the second 
derivative which meet the positive eigenvalue criterion men- 
tioned above. Third, we set the density derivatives r/ a i equal 
to zero. The last step is to read off the remainder of the cor- 
relation matrix for the velocities, much like what was done 
to derive Eq. 1)100 above. 
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